Bulk Density, Carbon, and OM
ds_cols <- names(depthseries)
# OM ~ DBD
if("fraction_organic_matter" %in% ds_cols & "dry_bulk_density" %in% ds_cols){
depthseries %>% drop_na(dry_bulk_density, fraction_organic_matter) %>%
ggplot(aes(x=fraction_organic_matter, dry_bulk_density, color = site_id)) +
geom_point(alpha = 0.6) +
xlab("Organic Matter (fraction)") +
ylab(expression(paste("Dry Bulk Density (g cm"^"-3",")", sep=""))) +
theme_bw() +
ggtitle("OM ~ DBD")
}
# C ~ DBD
if("fraction_carbon" %in% ds_cols & "dry_bulk_density" %in% ds_cols){
depthseries %>% drop_na(dry_bulk_density, fraction_carbon) %>%
ggplot(aes(x=fraction_carbon, dry_bulk_density, color = site_id)) +
geom_point(alpha = 0.6) +
xlab("Carbon (fraction)") +
ylab(expression(paste("Dry Bulk Density (g cm"^"-3",")", sep=""))) +
theme_bw() +
ggtitle("C ~ DBD")
}
# visual inspection to determine whether there are modeled values
if("fraction_carbon" %in% ds_cols & "fraction_organic_matter" %in% ds_cols){
depthseries %>% drop_na(fraction_organic_matter, fraction_carbon) %>%
ggplot(aes(x = fraction_organic_matter, y = fraction_carbon, color=site_id)) +
geom_point(alpha = 0.6) +
xlab("Organic Matter (fraction)") +
ylab("Organic Carbon (fraction)") +
theme_bw() +
ggtitle("Carbon ~ Organic Matter")
}
Depth Profiles
# Burden dataset problematic b/c there's only one 0-30 increment per core
if("fraction_organic_matter" %in% ds_cols){
depthseries %>% drop_na(fraction_organic_matter) %>%
ggplot(aes(x=depth_min, y=fraction_organic_matter)) +
geom_point(size = 2, pch=21, fill="white") +
geom_line() +
facet_wrap(~core_id, dir = "v") +
scale_x_reverse() +
ylab("Organic Matter (fraction)") +
xlab("Max Depth (cm)") +
coord_flip() +
theme_bw(base_size = 15) +
ggtitle("OM (LOI) Depth Profiles")
}
if("fraction_carbon" %in% ds_cols){
plot_c_profiles <- depthseries %>% drop_na(fraction_carbon) %>%
ggplot(aes(x=depth_min, y=fraction_carbon)) +
geom_line() +
geom_point(size = 2, pch=21, fill="white") +
facet_wrap(~core_id, dir = "v") +
scale_x_reverse() +
ylab("Carbon (fraction)") +
xlab("Max Depth (cm)") +
coord_flip() +
theme_bw(base_size = 15) +
ggtitle("C Depth Profiles")
}
This produces the plot with a special css class
if("fraction_carbon" %in% ds_cols){
plot_c_profiles
}
# there should be DBD, shouldn't need if statement
depthseries %>% drop_na(dry_bulk_density) %>%
ggplot(aes(x=depth_min, y=dry_bulk_density)) +
geom_line() +
geom_point(size = 2, pch=21, fill="white") +
facet_wrap(~core_id) +
ylab(expression(paste("Dry Bulk Density (g cm"^"-3",")", sep=""))) +
xlab("Max Depth (cm)") +
scale_x_reverse() +
coord_flip() +
theme_bw(base_size = 15) +
ggtitle("Bulk Density Profiles")
if(any(c("ra226_activity", "total_pb210_activity",
"excess_pb210_activity", "pb214_activity",
"bi214_activity") %in% ds_cols)){ # make this better
# total_pb210_activity_se may need to be calculated:
# pbsubset <- depthseries %>%
# drop_na(excess_pb210_activity, ra226_activity) %>%
# mutate(total_pb210_activity_se = excess_pb210_activity + ra226_activity)
# though it may already exist for some datasets
possible_ses <-c("ra226_activity_se", "total_pb210_activity_se",
"excess_pb210_activity_se", "pb214_activity_se",
"bi214_activity_se")
ses_there <- F
if (any(possible_ses%in%names(depthseries))) {
ses_there <- T
ses_in_dataset <- possible_ses[which(possible_ses%in%names(depthseries))]
}
possible_radioisotopes <-c("ra226_activity", "total_pb210_activity",
"excess_pb210_activity", "pb214_activity",
"bi214_activity")
radioisotoopes_in_dataset <- possible_radioisotopes[which(possible_radioisotopes%in%names(depthseries))]
if (ses_there) {
pbdata_se <- depthseries %>%
# mutate(ra226_activity_se = NA,
# total_pb210_activity_se = NA) %>%
select(core_id, depth_min,
ses_in_dataset)
names(pbdata_se) <- str_remove_all(names(pbdata_se), "_se")
pbdata_se <- pbdata_se %>%
gather(key="element", value="radioactivity_se", -core_id, -depth_min) %>%
drop_na(radioactivity_se)
}
pbdata <- depthseries %>%
select(core_id, depth_min, radioisotoopes_in_dataset) %>%
gather(key="element", value="radioactivity", -core_id, -depth_min) %>%
drop_na(radioactivity)
if (ses_there) {
pbdata <- pbdata %>%
left_join(pbdata_se)
}
pb_fig <- ggplot(pbdata, aes(x = depth_min, y = radioactivity, color = element)) +
geom_line() +
geom_point(size = 2, pch=21, fill="white") +
facet_wrap(~core_id) +
xlab("Max Depth (cm)") +
ylab(expression(paste("Radioactivity (dpm g"^"-1",")",sep=""))) +
scale_x_reverse() +
coord_flip() +
theme_bw(base_size = 15) +
theme(legend.position="bottom") +
ggtitle("Pb-210 Depth Profiles")
if (ses_there) {
pb_fig <- pb_fig + geom_segment(aes(x = depth_min, y = radioactivity - radioactivity_se,
xend = depth_min, yend = radioactivity + radioactivity_se))
}
(pb_fig)
}
if("age" %in% ds_cols){
ggplot(depthseries %>% drop_na(age), aes(x = depth_min, y = age)) +
geom_line() +
facet_wrap(~core_id) +
scale_x_reverse() +
ylab("Age") +
xlab("Max Depth (cm)") +
coord_flip() + theme_bw(base_size = 15) +
ggtitle("Age Depth Models") # plot error if present
}
if("cs137_activity" %in% ds_cols){
depthseries %>% drop_na(cs137_activity) %>%
ggplot(aes(x = depth_min, y = cs137_activity)) +
geom_line() +
geom_point(size = 2, pch=21, fill="white") +
facet_wrap(~core_id) +
scale_x_reverse() +
ylab("Cs-137 Activity") +
xlab("Max Depth (cm)") +
coord_flip() + theme_bw(base_size = 15) +
ggtitle("Cesium-137 Depth Profiles") # plot error if present
}
